Landcover Notebook¶
This notebook is being used to analyze the landcover in the area surrounding Lake Malawi. The work relies on data from the Sentinel-2 earth observation mission.
Dependencies¶
Here we'll import dependencies that can be used in other cells in this notebook.
import folium
import json
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import numpy as np
import pandas as pd
import pyproj
import rasterio
from rasterio.plot import show
from rasterio.enums import Resampling
from rasterio.mask import mask
from shapely.geometry import box, shape
from shapely.ops import transform
Helpers¶
Here we'll define some helper methods that can be used in other cells in this notebook.
classes = [1, 2, 4, 5, 7, 8, 9, 10, 11]
colors = [
"#419bdf", #1
"#397d49", #2
"#7a87c6", #4
"#e49635", #5
"#c4281b", #7
"#a59b8f", #8
"#a8ebff", #9
"#616161", #10
"#e3e2c3", #11
]
labels = [
"Water", #1
"Trees", #2
"Flooded", #4
"Crops", #5
"Built Area", #7
"Bare Ground", #8
"Snow/Ice", #9
"Clouds", #10
"Rangeland" #11
]
years = [2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024]
land_cover_colors = {
0: '#ffffff',
1: '#419bdf',
2: '#397d49',
4: '#7a87c6',
5: '#e49635',
7: '#c4281b',
8: '#a59b8f',
9: '#a8ebff',
10: '#616161',
11: '#e3e2c3',
}
land_cover_labels = {
1: 'Water',
2: 'Trees',
4: 'Flooded',
5: 'Crops',
7: 'Built Area',
8: 'Bare Ground',
9: 'Snow/Ice',
10: 'Clouds',
11: 'Rangeland',
}
class_label_map = {
"Class 1": "Water",
"Class 2": "Trees",
"Class 4": "Flooded",
"Class 5": "Crops",
"Class 7": "Built Area",
"Class 8": "Bare Ground",
"Class 9": "Snow/Ice",
"Class 10": "Clouds",
"Class 11": "Rangeland",
}
# File helpers
def input_path(filename):
return f"../data/input/{filename}"
def output_path(filename):
return f"../data/output/{filename}"
# Coordinate helpers
def get_crs(path):
with rasterio.open(path) as file:
return file.crs
def peek_coordinates(coordinates):
# Peek at the data
first_three_coordinates = coordinates[:3]
last_three_coordinates = coordinates[-3:]
print(f"{first_three_coordinates}...{last_three_coordinates}")
# Raster helpers
def create_mask(file, type, coords):
geometry = {
"type": type,
"coordinates": coords
}
return mask(file, [geometry], crop=True, nodata=file.nodata)
def get_stats_from_mask(masked_data, file):
valid_data = masked_data[masked_data != file.nodata]
unique_values, counts = np.unique(valid_data, return_counts=True)
stats = {
"classes": {},
"num_classes": len(unique_values),
"num_pixels": len(valid_data),
}
for value, count in zip(unique_values, counts):
percentage = (count / len(valid_data)) * 100
stats["classes"][f"Class {int(value)}"] = {
"pixels": count,
"percentage": percentage
}
return stats
def get_dataframe_from_mask(masked_data, file):
valid_data = masked_data[masked_data != file.nodata]
unique_values, counts = np.unique(valid_data, return_counts=True)
data = []
for value, num_pixels in zip(unique_values, counts):
percent_coverage = (num_pixels / len(valid_data)) * 100
data.append({
"Land Cover Class": int(value),
"Land Cover Label": land_cover_labels[int(value)],
"Number of Pixels": num_pixels,
"Percent Coverage": percent_coverage,
})
return pd.DataFrame(data)
Lake Malawi Coordinates¶
We need to get the coordinates of the perimiter of the Lake Malawi. For this, we'll use a geojson file obtained from PyGeoAPI.
with open(input_path("lake_malawi.json"), 'r') as f:
malawi_data = json.load(f)
malawi_coordinates = malawi_data['geometry']['coordinates'][0][0]
peek_coordinates(malawi_coordinates)
[35.2602047055542, -14.277474460510291]...[35.2602047055542, -14.277474460510291]
Map Lake Malawi¶
Lets take a peek at Lake Malawi.
map = folium.Map(location=[-11.6701, 34.6857], zoom_start=7)
folium.GeoJson(malawi_data).add_to(map)
map
Coordinate Reference System¶
In order to buffer the coordinates for Lake Malawi by 25km around the permiter and subsequently mask the geotiff files, we need to find out what coordinate reference system (CRS) we are working with. The CRS is embedded in the geotiff files.
crs = None
for year in years:
crs_from_file = get_crs(input_path(f"malawi_{year}.tif"))
if crs and crs != crs_from_file:
raise "File has unexpected CRS"
crs = crs_from_file
print(crs)
EPSG:32736
Buffered Lake Malawi Coordinates¶
Using the Lake Malawi coordinates and the known CRS, we will define a 25km buffer, defining new coordinates for the region of interest, which will be the lake and the 25km surrounding it.
malawi_geometry = shape(malawi_data['geometry'])
# Define coordinate transformation
# WGS84 (lat/lon) to a projected CRS so we can work accurately with area
transformer_to_utm = pyproj.Transformer.from_crs("EPSG:4326", crs, always_xy=True)
transformer_to_wgs84 = pyproj.Transformer.from_crs(crs, "EPSG:4326", always_xy=True)
# Transform to UTM (meters)
malawi_utm = transform(transformer_to_utm.transform, malawi_geometry)
# Buffer by 25km (25000 meters)
malawi_buffered = malawi_utm.buffer(25000)
malawi_buffered_coordinates = [list(malawi_buffered.exterior.coords)]
# Get the coordinates of the buffered perimter
malawi_buffered_wgs84 = transform(transformer_to_wgs84.transform, malawi_buffered)
malawi_buffered_wgs84_coordinates = [list(malawi_buffered_wgs84.exterior.coords)]
peek_coordinates(malawi_buffered_wgs84_coordinates)
[[(35.47707019386036, -14.19831045178232), (35.48353826652505, -14.217793727933875), (35.488171765266706, -14.237768211671032), (35.490932165245425, -14.258069838171183), (35.49179629546826, -14.278531816488693), (35.49075654464328, -14.298985995982939), (35.48782094029584, -14.319264245254097), (35.463971413975486, -14.44308510613051), (35.45876973377491, -14.464348288505134), (35.451481982267545, -14.485021816654664), (35.44217586254426, -14.504911914042387), (35.43093805487944, -14.523832086055018), (35.41787342399601, -14.541604873149732), (35.40310405270016, -14.55806352140359), (35.386768110284095, -14.57305355451169), (35.369018565791535, -14.586434232128523), (35.3500217578672, -14.598079880450072), (35.32995583444355, -14.607881082088564), (35.30900907693561, -14.615745713587659), (35.28737812489372, -14.621599820352737), (35.26526611818628, -14.625388320315738), (35.24288077472824, -14.627075529301862), (35.22043242252222, -14.626645502797926), (35.19813200532235, -14.624102190621302), (35.176189081556174, -14.619469402833264), (35.1548098362417, -14.612790587108915), (35.13419512551086, -14.60412841964682), (35.11453857300012, -14.593564213551373), (35.096024736798704, -14.581197150428274), (35.078827364864686, -14.567143342677516), (35.06310775584146, -14.551534735627946), (34.90538347947306, -14.378731176916634), (34.898021604340556, -14.389270579967956), (34.883933154733114, -14.40740268458517), (34.86804335055688, -14.424061405362687), (34.85051310109111, -14.439077558964824), (34.83152008799558, -14.452298594382686), (34.811256964777115, -14.463590153888948), (34.789929398386, -14.472837449718307), (34.76775397300152, -14.479946441869515), (34.744955977724, -14.484844804459346), (34.72176710130998, -14.487482670251419), (34.69842305824061, -14.487833145308242), (34.67516117128016, -14.485892588146077), (34.65221793624258, -14.481680650279786), (34.629826594930044, -14.475240077597276), (34.60821474213054, -14.46663627456686), (34.587601992164394, -14.455956635824025), (34.56819772975758, -14.44330965217374), (34.5501989690071, -14.4288238004478), (34.533788342907684, -14.412646228945727), (34.51913224435383, -14.39494125233344), (34.3594056450187, -14.180456038477432), (34.34777497928374, -14.163210656596632), (34.337807491986354, -14.14500102615036), (34.32958610857489, -14.125979342653716), (34.32317910417602, -14.106304555147533), (34.3186395453922, -14.086141032824782), (34.316004861194344, -14.065657187508757), (34.3152965460028, -14.045024063657962), (34.31902311088686, -13.747788769412766), (34.13780347169841, -13.517056452360542), (34.12510519979362, -13.49913765915147), (34.114232500643084, -13.4801041817296), (34.105286306772065, -13.46013344554673), (34.098349524532466, -13.439411573213372), (34.09348627654869, -13.418131643416071), (34.09074132189991, -13.396491886519334), (34.09013965875911, -13.374693833894641), (34.09859650490662, -12.985555710842416), (33.81881322478447, -12.291735467978812), (33.81177851879637, -12.271475577782951), (33.80672519734203, -12.250651926875786), (33.80369815094548, -12.229451413121623), (33.802724058739045, -12.208064273716827), (33.80381116623153, -12.186682375942446), (33.80694922900565, -12.16549749489286), (33.81210962221166, -12.144699593647493), (33.819245614251805, -12.12447512127258), (33.82829280161368, -12.105005343826205), (34.090527112586294, -11.60376106314986), (34.03441669001991, -10.510275908787255), (33.70605460586311, -9.90928658736391), (33.69648816532653, -9.889495160736555), (33.68890307589478, -9.868877014122889), (33.68337049736507, -9.847626575279122), (33.67994219993334, -9.825944198848143), (33.6786500903565, -9.804034274150649), (33.679505926485035, -9.78210329635905), (33.682501222261855, -9.76035791930653), (33.68760734347931, -9.739003008298166), (33.77678454047671, -9.432756716951495), (33.78413092874566, -9.411522725829501), (33.79356096826266, -9.39111675230756), (33.804981366812996, -9.371739878603599), (33.8182792722237, -9.353583010858616), (33.833323393872696, -9.336825003081884), (33.84996530309474, -9.32163090065757), (33.8680408994921, -9.308150320466), (33.887372028635205, -9.296515983320045), (33.907768235256974, -9.286842412917917), (33.929028634813825, -9.279224813880974), (33.95094388521654, -9.27373813969657), (33.973298239633536, -9.270436359535124), (33.9958716605439, -9.269351930976939), (34.01844197467468, -9.270495483685252), (34.04078704809925, -9.273855717017094), (34.062686960605284, -9.279399512492724), (34.083926158461026, -9.287072259966884), (34.10429556491837, -9.296798394282142), (34.123594628186645, -9.30848213715483), (34.14163328718938, -9.322008437066895), (34.158233836171064, -9.337244098031434), (34.68682462907354, -9.871693078503695), (34.70193302171426, -9.888530589016513), (34.71527679482443, -9.906775614784785), (34.72672320118735, -9.926246974316832), (34.73615826951539, -9.946751289617549), (34.74348794317016, -9.968084899809604), (34.74863902425353, -9.990035877318713), (34.75155991339243, -10.012386126940648), (34.830095664844436, -10.989942621399157), (35.11151579785281, -11.317093817062567), (35.12478248217734, -11.33405650794552), (35.13636288749836, -11.352182398464574), (35.14615274217831, -11.371308701129344), (35.15406379565459, -11.39126362372068), (35.16002461883554, -11.411867906149832), (35.163981255988865, -11.432936424924845), (35.16589772196135, -11.45427985108674), (35.16575633986412, -11.4757063469789), (35.163557915706996, -11.497023286832878), (35.15932174786272, -11.518038985899372), (34.927564431730474, -12.435298009560695), (35.07515351616887, -13.515539069363857), (35.10700689877398, -13.522599627525736), (35.1290484143779, -13.528620333738987), (35.150373648163445, -13.536744382406662), (35.1707736337812, -13.546892097196178), (35.19004841565827, -13.558963974687906), (35.20800899576726, -13.572841655077823), (35.224479176305806, -13.588389076307857), (35.239297280741525, -13.605453800579271), (35.2523177367293, -13.623868500574247), (35.263412505603945, -13.643452591195095), (35.27247234448313, -13.664013991240132), (35.47707019386036, -14.19831045178232)]]...[[(35.47707019386036, -14.19831045178232), (35.48353826652505, -14.217793727933875), (35.488171765266706, -14.237768211671032), (35.490932165245425, -14.258069838171183), (35.49179629546826, -14.278531816488693), (35.49075654464328, -14.298985995982939), (35.48782094029584, -14.319264245254097), (35.463971413975486, -14.44308510613051), (35.45876973377491, -14.464348288505134), (35.451481982267545, -14.485021816654664), (35.44217586254426, -14.504911914042387), (35.43093805487944, -14.523832086055018), (35.41787342399601, -14.541604873149732), (35.40310405270016, -14.55806352140359), (35.386768110284095, -14.57305355451169), (35.369018565791535, -14.586434232128523), (35.3500217578672, -14.598079880450072), (35.32995583444355, -14.607881082088564), (35.30900907693561, -14.615745713587659), (35.28737812489372, -14.621599820352737), (35.26526611818628, -14.625388320315738), (35.24288077472824, -14.627075529301862), (35.22043242252222, -14.626645502797926), (35.19813200532235, -14.624102190621302), (35.176189081556174, -14.619469402833264), (35.1548098362417, -14.612790587108915), (35.13419512551086, -14.60412841964682), (35.11453857300012, -14.593564213551373), (35.096024736798704, -14.581197150428274), (35.078827364864686, -14.567143342677516), (35.06310775584146, -14.551534735627946), (34.90538347947306, -14.378731176916634), (34.898021604340556, -14.389270579967956), (34.883933154733114, -14.40740268458517), (34.86804335055688, -14.424061405362687), (34.85051310109111, -14.439077558964824), (34.83152008799558, -14.452298594382686), (34.811256964777115, -14.463590153888948), (34.789929398386, -14.472837449718307), (34.76775397300152, -14.479946441869515), (34.744955977724, -14.484844804459346), (34.72176710130998, -14.487482670251419), (34.69842305824061, -14.487833145308242), (34.67516117128016, -14.485892588146077), (34.65221793624258, -14.481680650279786), (34.629826594930044, -14.475240077597276), (34.60821474213054, -14.46663627456686), (34.587601992164394, -14.455956635824025), (34.56819772975758, -14.44330965217374), (34.5501989690071, -14.4288238004478), (34.533788342907684, -14.412646228945727), (34.51913224435383, -14.39494125233344), (34.3594056450187, -14.180456038477432), (34.34777497928374, -14.163210656596632), (34.337807491986354, -14.14500102615036), (34.32958610857489, -14.125979342653716), (34.32317910417602, -14.106304555147533), (34.3186395453922, -14.086141032824782), (34.316004861194344, -14.065657187508757), (34.3152965460028, -14.045024063657962), (34.31902311088686, -13.747788769412766), (34.13780347169841, -13.517056452360542), (34.12510519979362, -13.49913765915147), (34.114232500643084, -13.4801041817296), (34.105286306772065, -13.46013344554673), (34.098349524532466, -13.439411573213372), (34.09348627654869, -13.418131643416071), (34.09074132189991, -13.396491886519334), (34.09013965875911, -13.374693833894641), (34.09859650490662, -12.985555710842416), (33.81881322478447, -12.291735467978812), (33.81177851879637, -12.271475577782951), (33.80672519734203, -12.250651926875786), (33.80369815094548, -12.229451413121623), (33.802724058739045, -12.208064273716827), (33.80381116623153, -12.186682375942446), (33.80694922900565, -12.16549749489286), (33.81210962221166, -12.144699593647493), (33.819245614251805, -12.12447512127258), (33.82829280161368, -12.105005343826205), (34.090527112586294, -11.60376106314986), (34.03441669001991, -10.510275908787255), (33.70605460586311, -9.90928658736391), (33.69648816532653, -9.889495160736555), (33.68890307589478, -9.868877014122889), (33.68337049736507, -9.847626575279122), (33.67994219993334, -9.825944198848143), (33.6786500903565, -9.804034274150649), (33.679505926485035, -9.78210329635905), (33.682501222261855, -9.76035791930653), (33.68760734347931, -9.739003008298166), (33.77678454047671, -9.432756716951495), (33.78413092874566, -9.411522725829501), (33.79356096826266, -9.39111675230756), (33.804981366812996, -9.371739878603599), (33.8182792722237, -9.353583010858616), (33.833323393872696, -9.336825003081884), (33.84996530309474, -9.32163090065757), (33.8680408994921, -9.308150320466), (33.887372028635205, -9.296515983320045), (33.907768235256974, -9.286842412917917), (33.929028634813825, -9.279224813880974), (33.95094388521654, -9.27373813969657), (33.973298239633536, -9.270436359535124), (33.9958716605439, -9.269351930976939), (34.01844197467468, -9.270495483685252), (34.04078704809925, -9.273855717017094), (34.062686960605284, -9.279399512492724), (34.083926158461026, -9.287072259966884), (34.10429556491837, -9.296798394282142), (34.123594628186645, -9.30848213715483), (34.14163328718938, -9.322008437066895), (34.158233836171064, -9.337244098031434), (34.68682462907354, -9.871693078503695), (34.70193302171426, -9.888530589016513), (34.71527679482443, -9.906775614784785), (34.72672320118735, -9.926246974316832), (34.73615826951539, -9.946751289617549), (34.74348794317016, -9.968084899809604), (34.74863902425353, -9.990035877318713), (34.75155991339243, -10.012386126940648), (34.830095664844436, -10.989942621399157), (35.11151579785281, -11.317093817062567), (35.12478248217734, -11.33405650794552), (35.13636288749836, -11.352182398464574), (35.14615274217831, -11.371308701129344), (35.15406379565459, -11.39126362372068), (35.16002461883554, -11.411867906149832), (35.163981255988865, -11.432936424924845), (35.16589772196135, -11.45427985108674), (35.16575633986412, -11.4757063469789), (35.163557915706996, -11.497023286832878), (35.15932174786272, -11.518038985899372), (34.927564431730474, -12.435298009560695), (35.07515351616887, -13.515539069363857), (35.10700689877398, -13.522599627525736), (35.1290484143779, -13.528620333738987), (35.150373648163445, -13.536744382406662), (35.1707736337812, -13.546892097196178), (35.19004841565827, -13.558963974687906), (35.20800899576726, -13.572841655077823), (35.224479176305806, -13.588389076307857), (35.239297280741525, -13.605453800579271), (35.2523177367293, -13.623868500574247), (35.263412505603945, -13.643452591195095), (35.27247234448313, -13.664013991240132), (35.47707019386036, -14.19831045178232)]]
Save Buffered Lake Malawi Coordinates¶
We will save a new geojson file with the buffered coordinates.
expanded_geojson = {
"type": "Feature",
"properties": {
"id": 10,
"scalerank": 0,
"name": "Lake Malawi",
"name_alt": "Lake Nyasa",
"admin": "admin-0",
"featureclass": "Lake"
},
"geometry": {
"type": malawi_buffered_wgs84.geom_type,
"coordinates": malawi_buffered_wgs84_coordinates
},
"id": 10
}
with open(output_path("lake_malawi_expanded_25km.json"), "w") as dst:
json.dump(expanded_geojson, dst, indent=2)
Map Lake Malawi and Surrounding 25KM¶
Map the expanded boundary of Lake Malawi.
folium.GeoJson(expanded_geojson).add_to(map)
map
Visualize Landcover for Greater Region¶
Let's visualize the landcover for the greater region using the original geotiff files.
for year in years:
max_class_value = max(land_cover_colors.keys())
cmap_colors = [land_cover_colors.get(i, '#ffffff') for i in range(max_class_value + 1)]
custom_cmap = ListedColormap(cmap_colors)
with rasterio.open(input_path(f"malawi_{year}.tif")) as file:
downsampled_data = file.read(
out_shape=(
file.count,
int(file.height * 0.1),
int(file.width * 0.1)
),
resampling=Resampling.bilinear
)
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
show(downsampled_data, ax=ax, transform=file.transform, title = f"Landcover in {year}", cmap=custom_cmap)
legend_labels = [f'{val}: {name}' for val, name in land_cover_labels.items()]
plt.legend(handles=[plt.Rectangle((0,0),1,1, color=c) for c in cmap_colors if c != '#ffffff'], labels=legend_labels, loc='lower left', bbox_to_anchor=(1.05, 0))
plt.axis('off')
plt.show()
Analyze Land Cover within Buffered Lake Malawi Coordinates¶
The geotiff files contain the land cover data for a large region surrounding Lake Malawi. We'll mask the files to get the cover for the immediate area surrounding the lake.
for year in years:
with rasterio.open(input_path(f"malawi_{year}.tif")) as file:
masked_data, masked_transform = create_mask(
file,
malawi_buffered.geom_type,
malawi_buffered_coordinates
)
bounding_box = box(*file.bounds)
if bounding_box.intersects(malawi_buffered):
print(f"\n{year}")
stats = get_stats_from_mask(masked_data, file)
print(f"Number of pixels: {stats["num_pixels"]}")
print(f"Number of land cover classes: {stats["num_classes"]}")
df = get_dataframe_from_mask(masked_data, file)
print(df)
else:
print("No overlap detected. Mask is not correct")
output_profile = file.profile.copy()
output_profile.update({
'height': masked_data.shape[1],
'width': masked_data.shape[2],
'transform': masked_transform
})
with rasterio.open(output_path(f"malawi_{year}_masked.tif"), 'w', **output_profile) as dst:
dst.write(masked_data)
2017 Number of pixels: 594349727 Number of land cover classes: 8 Land Cover Class Land Cover Label Number of Pixels Percent Coverage 0 1 Water 296062405 49.812828 1 2 Trees 94756092 15.942817 2 4 Flooded 279097 0.046958 3 5 Crops 23663623 3.981431 4 7 Built Area 8954681 1.506635 5 8 Bare Ground 163789 0.027558 6 10 Clouds 5742 0.000966 7 11 Rangeland 170464298 28.680807 2018 Number of pixels: 594349727 Number of land cover classes: 8 Land Cover Class Land Cover Label Number of Pixels Percent Coverage 0 1 Water 296140745 49.826008 1 2 Trees 110617108 18.611451 2 4 Flooded 544220 0.091566 3 5 Crops 22283444 3.749214 4 7 Built Area 10205910 1.717156 5 8 Bare Ground 96417 0.016222 6 10 Clouds 84383 0.014198 7 11 Rangeland 154377500 25.974185 2019 Number of pixels: 594349727 Number of land cover classes: 8 Land Cover Class Land Cover Label Number of Pixels Percent Coverage 0 1 Water 296188895 49.834110 1 2 Trees 110611997 18.610591 2 4 Flooded 458157 0.077085 3 5 Crops 23421564 3.940704 4 7 Built Area 10550809 1.775185 5 8 Bare Ground 83680 0.014079 6 10 Clouds 25875 0.004353 7 11 Rangeland 153008750 25.743892 2020 Number of pixels: 594349727 Number of land cover classes: 8 Land Cover Class Land Cover Label Number of Pixels Percent Coverage 0 1 Water 296398142 49.869316 1 2 Trees 100409978 16.894090 2 4 Flooded 816341 0.137350 3 5 Crops 25856763 4.350429 4 7 Built Area 10615203 1.786020 5 8 Bare Ground 57904 0.009742 6 10 Clouds 6062 0.001020 7 11 Rangeland 160189334 26.952033 2021 Number of pixels: 594349727 Number of land cover classes: 8 Land Cover Class Land Cover Label Number of Pixels Percent Coverage 0 1 Water 296437041 49.875861 1 2 Trees 86120127 14.489807 2 4 Flooded 916144 0.154142 3 5 Crops 26686570 4.490045 4 7 Built Area 11680205 1.965207 5 8 Bare Ground 46410 0.007809 6 10 Clouds 16203 0.002726 7 11 Rangeland 172447027 29.014403 2022 Number of pixels: 594349727 Number of land cover classes: 8 Land Cover Class Land Cover Label Number of Pixels Percent Coverage 0 1 Water 296576042 49.899248 1 2 Trees 95916954 16.138134 2 4 Flooded 1297987 0.218388 3 5 Crops 30356672 5.107544 4 7 Built Area 12271820 2.064747 5 8 Bare Ground 38917 0.006548 6 10 Clouds 22781 0.003833 7 11 Rangeland 157868554 26.561559 2023 Number of pixels: 594349727 Number of land cover classes: 8 Land Cover Class Land Cover Label Number of Pixels Percent Coverage 0 1 Water 296870362 49.948767 1 2 Trees 94974139 15.979504 2 4 Flooded 2233678 0.375819 3 5 Crops 27576826 4.639832 4 7 Built Area 12267058 2.063946 5 8 Bare Ground 28997 0.004879 6 10 Clouds 24131 0.004060 7 11 Rangeland 160374536 26.983193 2024 Number of pixels: 594349727 Number of land cover classes: 8 Land Cover Class Land Cover Label Number of Pixels Percent Coverage 0 1 Water 297201661 50.004509 1 2 Trees 103221610 17.367150 2 4 Flooded 2786235 0.468787 3 5 Crops 31337596 5.272585 4 7 Built Area 12553491 2.112139 5 8 Bare Ground 35388 0.005954 6 10 Clouds 4004 0.000674 7 11 Rangeland 147209742 24.768202
Visualize Landcover for Lake Malawi and 25KM Surroudning It¶
Let's visualize the landcover for the buffered Lake Malawi area.
for year in years:
max_class_value = max(land_cover_colors.keys())
cmap_colors = [land_cover_colors.get(i, '#ffffff') for i in range(max_class_value + 1)]
custom_cmap = ListedColormap(cmap_colors)
with rasterio.open(output_path(f"malawi_{year}_masked.tif")) as file:
downsampled_data = file.read(
out_shape=(
file.count,
int(file.height * 0.1),
int(file.width * 0.1)
),
resampling=Resampling.bilinear
)
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
image = show(downsampled_data, ax=ax, transform=file.transform, title = f"Landcover in {year} (25km around Lake Malawi)", cmap=custom_cmap)
legend_labels = [f'{val}: {name}' for val, name in land_cover_labels.items()]
plt.legend(handles=[plt.Rectangle((0,0),1,1, color=c) for c in cmap_colors if c != '#ffffff'], labels=legend_labels, loc='lower left', bbox_to_anchor=(1.05, 0))
plt.axis('off')
plt.show()
Visualize Landcover Changes¶
See the changes in landcover from 2017-2024
frames = []
for year in years:
with rasterio.open(input_path(f"malawi_{year}.tif")) as file:
masked_data, masked_transform = create_mask(
file,
malawi_buffered.geom_type,
malawi_buffered_coordinates
)
bounding_box = box(*file.bounds)
if bounding_box.intersects(malawi_buffered):
df = get_dataframe_from_mask(masked_data, file)
frames.append(df)
for i, year in enumerate(years):
frames[i]['Year'] = year
combined_df = pd.concat(frames, ignore_index=True)
print(combined_df)
pivot_df = combined_df.pivot(index='Year', columns='Land Cover Label', values='Percent Coverage')
fig, ax = plt.subplots(figsize=(12, 6))
for label in pivot_df.columns:
class_val = [k for k, v in land_cover_labels.items() if v == label][0]
color = land_cover_colors.get(class_val, '#000000')
ax.plot(pivot_df.index, pivot_df[label], marker='o', label=label, color=color, linewidth=2)
ax.set_xlabel('Year', fontsize=12)
ax.set_ylabel('Percent Coverage (%)', fontsize=12)
ax.set_title('Land Cover Change Over Time (2017-2024)', fontsize=14)
ax.legend(loc='best')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
Land Cover Class Land Cover Label Number of Pixels Percent Coverage \
0 1 Water 296062405 49.812828
1 2 Trees 94756092 15.942817
2 4 Flooded 279097 0.046958
3 5 Crops 23663623 3.981431
4 7 Built Area 8954681 1.506635
.. ... ... ... ...
59 5 Crops 31337596 5.272585
60 7 Built Area 12553491 2.112139
61 8 Bare Ground 35388 0.005954
62 10 Clouds 4004 0.000674
63 11 Rangeland 147209742 24.768202
Year
0 2017
1 2017
2 2017
3 2017
4 2017
.. ...
59 2024
60 2024
61 2024
62 2024
63 2024
[64 rows x 5 columns]
Ideas for Next Up¶
- Get evapotranspiration data from https://planetarycomputer.microsoft.com/explore?c=-105.1007%2C40.1705&z=7.00&v=2&d=modis-16A3GF-061&m=Most+recent&r=Total+evapotranspiration+yearly+500m+%28kg%2Fm²%2Fyear%29&s=false%3A%3A100%3A%3Atrue&sr=desc&ae=0
- Clamp data to the 25km buffer for Lake Malawi
- Visualize data
- Visualize change over time
- Compare change over time to changes in landcover
- Collect addition measurements for 25km buffer for Lake Malawi
- ML using random decision forest using landcover + additional measurements to predict evoptransipration
- Assess accuracy of prediction
- Try to find most important factors